#f Author: Ivan Canay | Northwestern University | Nov 2019
# Copyright (c) 2019 Northwestern University. All rights reserved.
#-------------------------------------------------------------------
# X-plain: Computes the simulations for Wild Bootstrap paper 
#------------------------------------------------------------------
rm(list=ls(all=TRUE))               # Clear all variables
source("Wild.Functions.R")
library(xtable)

#-------------------------------------------------------------------
Wild.Break <-function(myseed,design,q,n,alpha=0.05,MC=5000,beta=1){
#----------------------------------------------------------------
#IMPUTS: myseed: seed for the simulations 
#		 design: 3-4 
#        q: number of clusters
#        fixed = 1 if the Z are fixed across simulations
#        beta: true value is equal to 1 under the null and 0 under alt.  
#        MC: Number of monte carlo simulations
# Other parameters

null.b = 1; # value tested under H_0
k= 3; # number of covariates

# GENERATE DATA and SET PARAMETERS
# - - - - - - - - - - - - - - - - -
dir.create("Tables")
nametext <-paste("Tables/Wild.Break",q,n,"design",design,myseed,"txt",sep=".");
dir.create("Latex")
nameLaTeX<-paste("Latex/Wild.Break",q,n,"design",design,myseed,"tex",sep=".");
dir.create("RData")
Rname 	 <-paste("RData/Wild.Break",q,n,"design",design,myseed,"Rdata",sep=".");
set.seed(myseed);

# Set all rejection probabilities to zero
# save test statistics and cv
rej.NS <-matrix(0,MC,4); rej.S <- rej.NS; rej.C <- rej.NS; 
rej.EN <- rej.NS; rej.ES <- rej.NS; rej.EC <- rej.NS;
T.NS.simple<-matrix(0,MC,1); T.NS.FE<-T.NS.simple; T.S.simple<-T.NS.FE; T.S.FE<-T.NS.FE;
cv.NS.rade.simple<-matrix(0,MC,1);
cv.NS.mame.simple<-cv.NS.rade.simple; 
cv.NS.rade.FE<-cv.NS.rade.simple; cv.NS.mame.FE<-cv.NS.rade.simple;
cv.S.rade.simple<-cv.NS.rade.simple;
cv.S.mame.simple<-cv.NS.rade.simple; 
cv.S.rade.FE<-cv.NS.rade.simple; cv.S.mame.FE<-cv.NS.rade.simple;
cv.C.rade.simple<-cv.NS.rade.simple;
cv.C.mame.simple<-cv.NS.rade.simple; 
cv.C.rade.FE<-cv.NS.rade.simple; cv.C.mame.FE<-cv.NS.rade.simple;

time1 <- proc.time()
time3 <- proc.time()

core = gen.corr(k,q);

if (design==4){
	k=2;
	core = gen.corr(k,q);
	C1=matrix(c(10,0.8,0.8,1),2,2);
	C2=diag(2);
	core[[1]] = C1;
	core[[2]] = C2;
}

for (j in 1:MC){
	data = gen.data.k(q,n,core,beta);
	clu  = data$cluster;
	N    = length(clu);

	# TEST STATISTICS - two regressions
	simple.reg 	<- lm(Y~ .-cluster-unit,data);
	simple.beta <- coef(simple.reg)[2];
	FE.reg    	<- lm(Y~ .-cluster-unit+factor(cluster),data);
	FE.beta 	<- coef(FE.reg)[2];

	simple.covariates  = gen.covariates(data,0);
	FE.covariates      = gen.covariates(data,1);

	simple.cce  = sqrt(N*cce.stata(simple.covariates,simple.reg,clu,"stata")[2,2]);
	FE.cce  	= sqrt(N*cce.stata(FE.covariates,FE.reg,clu,"stata")[2,2]);

	T.NS.simple[j] <- sqrt(N)*(simple.beta-null.b);
	T.NS.FE[j]     <- sqrt(N)*(FE.beta-null.b);

	T.S.simple[j]  <- sqrt(N)*(simple.beta-null.b)/simple.cce;
	T.S.FE[j] 	   <- sqrt(N)*(FE.beta-null.b)/FE.cce;

	# RUN ALL THE WILD BOOTSTRAPS
	bo.rade.simple = Wild.Boot.k(data,null.b,alpha,0);
	bo.rade.FE 	   = Wild.Boot.k(data,null.b,alpha,1);

	# CRITICAL VALUES FOR THE ABSOLUTE VALUE OF TEST STATISTIC
	# Non-studentized
	cv.NS.rade.simple[j] <- bo.rade.simple$cv.NS;
	cv.NS.rade.FE[j] 	 <- bo.rade.FE$cv.NS;

	# Studentized
	cv.S.rade.simple[j]  <- bo.rade.simple$cv.S;
	cv.S.rade.FE[j]  	 <- bo.rade.FE$cv.S;

	# Studentized corrected
	cv.C.rade.simple[j]  <- bo.rade.simple$cv.Sco;
	cv.C.rade.FE[j]  	 <- bo.rade.FE$cv.Sco;

	# EQUI-TAIL CRITICAL VALUES
	# Non-studentized
	ET.NS.rade.Si <- c(bo.rade.simple$cL.NS,bo.rade.simple$cH.NS);
	ET.NS.rade.FE <- c(bo.rade.FE$cL.NS,bo.rade.FE$cH.NS);

	# Studentized
	ET.S.rade.Si  <- c(bo.rade.simple$cL.S,bo.rade.simple$cH.S);
	ET.S.rade.FE  <- c(bo.rade.FE$cL.S,bo.rade.FE$cH.S);

	# Studentized Corrected
	ET.C.rade.Si  <- c(bo.rade.simple$cL.Sco,bo.rade.simple$cH.Sco);
	ET.C.rade.FE  <- c(bo.rade.FE$cL.Sco,bo.rade.FE$cH.Sco);

	# Compute rejection probabilities of all tests - ABSOLUTE VALUE
	rej.NS.rade.simple <-(abs(T.NS.simple[j])>cv.NS.rade.simple[j]);
	rej.NS.rade.FE     <-(abs(T.NS.FE[j])>cv.NS.rade.FE[j]);
	rej.S.rade.simple <-(abs(T.S.simple[j])>cv.S.rade.simple[j]);
	rej.S.rade.FE     <-(abs(T.S.FE[j])>cv.S.rade.FE[j]);
	
	rej.C.rade.simple <-(abs(T.S.simple[j])>=cv.C.rade.simple[j]); # note rejection region is "">=""
	rej.C.rade.FE     <-(abs(T.S.FE[j])>=cv.C.rade.FE[j]);

	# Compute rejection probabilities of all tests - EQUI TAILED
	rej.EN.rade.simple<- (T.NS.simple[j]<ET.NS.rade.Si[1]) + (T.NS.simple[j]>ET.NS.rade.Si[2]);
	rej.EN.rade.FE 	  <- (T.NS.FE[j]<ET.NS.rade.FE[1]) + (T.NS.FE[j]>ET.NS.rade.FE[2]);

	rej.ES.rade.simple<- (T.S.simple[j]<ET.S.rade.Si[1]) + (T.S.simple[j]>ET.S.rade.Si[2]);
	rej.ES.rade.FE    <- (T.S.FE[j]<ET.S.rade.FE[1]) + (T.S.FE[j]>ET.S.rade.FE[2]);


	rej.EC.rade.simple<- (T.S.simple[j]<=ET.C.rade.Si[1]) + (T.S.simple[j]>=ET.C.rade.Si[2]);
	rej.EC.rade.FE    <- (T.S.FE[j]<=ET.C.rade.FE[1]) + (T.S.FE[j]>=ET.C.rade.FE[2]);

	# PUT ALL REJECTIONS TOGETHER
	rej.NS[j,] <- c(rej.NS.rade.simple,rej.NS.rade.FE,0,0);
	rej.S[j,]  <- c(rej.S.rade.simple,rej.S.rade.FE,0,0);
	rej.C[j,]  <- c(rej.C.rade.simple,rej.C.rade.FE,0,0);
	rej.EN[j,] <- c(rej.EN.rade.simple,rej.EN.rade.FE,0,0);
	rej.ES[j,] <- c(rej.ES.rade.simple,rej.ES.rade.FE,0,0);
	rej.EC[j,] <- c(rej.EC.rade.simple,rej.EC.rade.FE,0,0);

	# Report Iterations on Screen for monitoring
	if (floor(10*j/MC)==10*j/MC){
    cat("% completed = ", j/MC*100, "\n")
	  time4 <- proc.time();
    cat("Time Last round (min) = ", (time4[3]-time3[3])/60, "\n")
    time3 <-proc.time()
	} 
}
time2 <- proc.time();
cat(sprintf("Elapsed Time (min): %1.3f",(time2[3]-time1[3])/60), "\n");

## Display Information
basics_a <- c(q,n,alpha,MC);
names(basics_a)<- c("   q","  n", "alpha", "MC simu");

## Null rejection probabilities
NonStudentized <- colMeans(rej.NS)*100;
Studentized    <- colMeans(rej.S)*100;
Corrected      <- colMeans(rej.C)*100;
Equi.NonStu    <- colMeans(rej.EN)*100;
Equi.Stu       <- colMeans(rej.ES)*100;
Equi.Corrected <- colMeans(rej.EC)*100;

names(NonStudentized)  <-c("                Simple Rade |", "  FE Rade  |", " Simple Mammen  |", "  FE Mammen |");
results <-rbind(NonStudentized,Studentized,Corrected,Equi.NonStu,Equi.Stu,Equi.Corrected);

# Save LaTeX table
X.results<-cbind(c("NonStudentized","Studentized","Corrected","Equi NonStu","Equi Stu","Equi Corrected"),results)
# VERY IMPORTANT!!! row names cannot be duplicates, so make sure there are different "spaces" in each 
# name that is an empty name. 
row.names(X.results)<-c(design," ","  ","   ","    ","");
print(X.results)
X.results<-xtable(X.results,digits=2);
print(X.results, type="latex", file=nameLaTeX);

print(results)

# # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Saving information to a text file
out <- file(nametext,"wt");   # Choose name of the file here
cat("---------------------------------------------------------------------------------------------------------------------", "\n",file=out);
cat(date(), "\n",file=out)
cat("Seed of the round = ", myseed , "\n",file=out)
cat("Design used       = ", design , "\n",file=out)
cat("Number of covariates = ", k , "\n",file=out)
cat("----------------------------------------------- Parameters ----------------------------------------------------------", "\n",file=out);
cat(names(basics_a),sep="   ","\n",file=out);
cat(formatC(basics_a,digits=2,format="f"),sep="   ", "\n",file=out);
cat("\n",file=out)
cat("---------------------------------------- REJECTION PROBABILITIES ----------------------------------------------------", "\n",file=out);
cat(names(NonStudentized),sep=" ","\n",file=out);
write.table(format(results, drop0trailing=FALSE),file=out,quote=F,sep='\t\t',row.names=T,col.names=F)
cat("\n",file=out)
cat("---------------------------------------------------------------------------------------------------------------------", "\n",file=out);
time2 <- proc.time();
cat(sprintf("Elapsed Time (min): %1.3f",(time2[3]-time1[3])/60), "\n",file=out);
cat("---------------------------------------------------------------------------------------------------------------------", "\n",file=out);
cat(" Simple: means no fixed effects. FE: means regression with fixed effects", "\n",file=out);
cat(" Equitail use the alpha/2 and 1-alpha/2 quantiles and rejects if T<q(alpha/2) or T>q(1-alpha/2)", "\n",file=out);
cat(" Design 1: Baseline design based on Jamie and Max ", "\n",file=out);
cat("---------------------------------------------------------------------------------------------------------------------", "\n",file=out);
close(out);

save(list = ls(all = TRUE), file = Rname)
}